rm(list=ls())
setwd("")
library(haven)


load("grid_data.Rdata")

##HARBOR VARIABLE 

uniquesetVdem$natharbdist_km <- uniquesetVdem$two_port_all_pred1/1000

summary(uniquesetVdem$natharbdist_km)

uniquesetVdem$natharbdist_km <- uniquesetVdem$natharbdist/100

summary(uniquesetVdem$natharbdist_km)

grid <- uniquesetVdem
rm(uniquesetvdem)

#new polyarchy
library(haven)
pol <- read_dta("Demo_book_country_3.dta")

names(pol)
pol <- subset(pol, select=c(country_id, year, v2x_polyarchy_imp))

#merging 
grid <- merge(grid, pol, by=c("country_id", "year"), all.x=T)
#################################################################################
############             INTRODUCING THE META GRID CELLS  ##########
#################################################################################
library(sp)
library(rgdal)

#4x2 decimal degrees 
grid1 <- readOGR("Making the grid-cell dataset/largegrids/4x2dd.shp", layer="4x2dd")
proj4string(grid1)
proj4string(e_priogrid)
e1 <- over(e_priogrid, grid1)
e1 <- cbind(e_priogrid@data$gid,e1)
colnames(e1) <- c("gid", "4x2dd")

#6x3 decimal degrees
grid2 <- readOGR("Making the grid-cell dataset/largegrids/6x3dd.shp", layer="6x3dd")
e2 <- over(e_priogrid, grid2)
e1 <- cbind(e1,e2)

#8x4 decimal degrees
grid3 <- readOGR("Making the grid-cell dataset/largegrids/8x4dd.shp", layer="8x4dd")
e3 <- over(e_priogrid, grid3)
e1 <- cbind(e1,e3)

#12x6 decimal degrees
grid4 <- readOGR("Making the grid-cell dataset/largegrids/12x6dd.shp", layer="12x6dd")
e4 <- over(e_priogrid, grid4)
e1 <- cbind(e1,e4)


#18x9 decimal degrees
grid5 <- readOGR("Making the grid-cell dataset/largegrids/18x9dd.shp", layer="18x9dd")
e5 <- over(e_priogrid, grid5)
e1 <- cbind(e1,e5)


#24x12 decimal degrees
grid6 <- readOGR("Making the grid-cell dataset/largegrids/24x12dd.shp", layer="24x12dd")
e6 <- over(e_priogrid, grid6)
e1 <- cbind(e1,e6)


#36x18dd decimal degrees
grid7 <- readOGR("Making the grid-cell dataset/largegrids/36x18dd.shp", layer="36x18dd")
e7 <- over(e_priogrid, grid7)
e1 <- cbind(e1,e7)


#72x36 decimal degrees
grid8 <- readOGR("Making the grid-cell dataset/largegrids/72x36dd.shp", layer="72x36dd")
e8 <- over(e_priogrid, grid8)
e1 <- cbind(e1,e8)

colnames(e1) <- c("gid", "4x2dd", "6x3dd", "8x4dd", "12x6dd", "18x9dd", "24x12dd", "36x18dd","72x36dd")

rm(e2,e3,e4,e5,e6,e7,e8,
   grid1,grid2,grid3,grid4,grid5,grid6,grid7,grid8,grid9, loyds)

#merge 
grid <- merge(grid, e1, by=c("gid"), all.x=T)



###AGGREGATE TO META_GRIDS
g1 <- subset(grid, select=c(v2x_polyarchy_imp, natharbdist_km, equator_distance
,prec_gpcc,rivdist , irrig_sum , temp , crop , pasture,
  forest ,grass ,shrub ,savanna , barren , water  , mountains, year))

g1 <- cbind(g1, grid[,117:124])
names(g1)
data <- g1

d1 <- aggregate(data, by=list(data[,18], data$year), FUN=mean, na.rm=T)
d2 <- aggregate(data, by=list(data[,19], data$year), FUN=mean, na.rm=T)
d3 <- aggregate(data, by=list(data[,20], data$year), FUN=mean, na.rm=T)
d4 <- aggregate(data, by=list(data[,21], data$year), FUN=mean, na.rm=T)
d5 <- aggregate(data, by=list(data[,22], data$year), FUN=mean, na.rm=T)
d6 <- aggregate(data, by=list(data[,23], data$year), FUN=mean, na.rm=T)
d7 <- aggregate(data, by=list(data[,24], data$year), FUN=mean, na.rm=T)
d8 <- aggregate(data, by=list(data[,25], data$year), FUN=mean, na.rm=T)


#regressions
#m1
summary(m0 <- lm(v2x_polyarchy_imp ~ natharbdist_km + as.factor(year)-1, data=data))
summary(m0)
#m1
summary(m1 <- lm(v2x_polyarchy_imp ~ natharbdist_km + as.factor(year)-1, data=d1))
#m2
summary(m2 <- lm(v2x_polyarchy_imp ~ natharbdist_km + as.factor(year)-1 ,data=d2))
#m3
summary(m3 <- lm(v2x_polyarchy_imp ~ natharbdist_km + as.factor(year)-1, data=d3))
#m4
summary(m4 <- lm(v2x_polyarchy_imp ~ natharbdist_km + as.factor(year)-1, data=d4))

#m5
summary(m5 <- lm(v2x_polyarchy_imp ~ natharbdist_km + as.factor(year)-1, data=d5))
#m6
summary(m6 <- lm(v2x_polyarchy_imp ~ natharbdist_km + as.factor(year)-1, data=d6))
#m7
summary(m7 <- lm(v2x_polyarchy_imp ~ natharbdist_km + as.factor(year)-1, data=d7))

#m8
summary(m8 <- lm(v2x_polyarchy_imp ~ natharbdist_km + as.factor(year)-1, data=d8))


coefs <- rbind(  coef(m2)[1:1],
                coef(m3)[1:1], coef(m4)[1:1], coef(m5)[1:1], coef(m6)[1:1], 
                coef(m7)[1:1], coef(m8)[1:1], coef(m0)[1:1])

ses<- rbind( coef(summary(m2))[1, "Std. Error"],
            coef(summary(m3))[1, "Std. Error"], coef(summary(m4))[1, "Std. Error"], 
            coef(summary(m5))[1, "Std. Error"],coef(summary(m6))[1, "Std. Error"],
            coef(summary(m7))[1, "Std. Error"], coef(summary(m8))[1, "Std. Error"],
            coef(summary(m0))[1, "Std. Error"],coef(summary(m0))[1, "Std. Error"])


#plot of coefficients

dds <- c( 6, 12, 
         18, 24, 36 , 72, 180)
dds <- (360/dds)*110
ns <- c(  3888, 5703, 11281, 21698, 32212, 61255, 183342)
dds <- paste(dds,"x", dds, "Km", "(N=",ns,")", sep="")
dds <- c(dds, "50x50km(N=10885183)")


tiff("metagrids.tiff", res=200)
par(mai=c(1,2.4,1,1))
plot(NA, xlim = c(-7, 0), ylim = c(0, 8), xlab = "Coefficient (Harbor distance)", ylab = "", yaxt = "n")
# We can add a title:
title("Estimated coefficients for different grid-cell-sizes")
# We'll add a y-axis labelling our variables:


axis(2, 1:8, dds, las = 2)
# We'll add a vertical line for zero:
abline(v = 0, col = "gray")
# Then we'll draw our slopes as points (`pch` tells us what type of point):
points(coefs[,1], 1:8, pch = 23, col = "grey", bg = "grey")

# Then we'll add thick line segments for each 1 SE:

for(i in 1:8){
  segments((coefs[i,1] - qnorm(0.975) * ses[i,1]), i, (coefs[i,1] + qnorm(0.975) *ses[i,1]), i, col = "black", lwd = 2)
}
dev.off()
